arXiv:1506.03288v3 [cond-mat.stat-mech] 21 Oct 2015 


Protocol-independent granular temperature supported by numerical simulations 

Volker BeckeiQ and Klaus KassneiQ 
Institute for Theoretical Physics, 

Otto von Guericke University Magdeburg 
Postfach 4130 , 39106 Magdeburg, Germany 
(Dated: August 27, 2015) 

A possible approach to the statistical description of granular assemblies starts from Edwards’ 
assumption that all blocked states occupying the same volume are equally probable (S.F. Ed¬ 
wards, R. Oakeshott, Physica A 157, 1080 (1989)). We performed computer simulations using 
two-dimensional polygonal particles excited periodically according to two different protocols: ex¬ 
citation by pulses of “negative gravity” and excitation by “rotating gravity”. The first protocol 
exhibits a non-monotonous dependency of the mean volume fraction on the pulse strength. The 
overlapping histogram method is used in order to test whether or not the volume distribution is 
described by a Boltzmann-like distribution, and to calculate the inverse compactivity as well as the 
logarithm of the partition sum. We find that the mean volume is a unique function of the mea¬ 
sured granular temperature, independently of the protocol and of the branch in and that all 

determined quantities are in agreement with Edwards’ theory. 

PACS numbers: 45.70.-n, 05.20.Gg 


I. INTRODUCTION 

Granular materials are typically composed of thou¬ 
sands to millions of individual particles (or more). Think 
for example of the cereals for breakfast or of the sand on 
the shore. These large numbers suggest that statistical 
methods may be applicable and constitute a powerful 
tool in developing a better theoretical understanding of 
these kinds of materials. However, contrary to the situa¬ 
tion in gases or fluids, where the permitted phase space is 
explored continuously due to chaotic molecular motion, 
thermal fluctuations are negligible for granular materi¬ 
als. Furthermore, the particle dynamics are dissipative. 
Therefore, it is not possible to carry standard statistical 
mechanics over to granular assemblies. In particular, the 
classical Boltzmann distribution, where the probability 
of a state is inversely proportional to the exponential of 
its energy (measured in units of ksT), will not apply to 
these systems. 

Edwards and Oakeshott [HE] proposed that concepts 
from classical statistical mechanics are applicable, if one 
assumes that the volume of a static, stable granulate 
plays the same role in granular statistics as the energy of 
a microstate in classical statistics. This means that, anal¬ 
ogous to the classical microcanonical ensemble, where all 
states with the same energy are equally probable, all me¬ 
chanically stable configurations of the granular assem¬ 
bly that occupy the same volume occur with the same 
probability. The entropy of this granular microcanonical 
ensemble is proportional to the logarithm of the num¬ 
ber of blocked states with a certain volume. By analogy 
with classical statistics, one can define a temperature¬ 
like variable, called compactivity, as the inverse of the 
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derivative of the Edwards entropy with respect to the 
volume. Later, it turned out that for a full description 
of a granular system, beyond the volume ensemble also 
an ensemble for the different stress states must be intro¬ 
duced mg and that the volume and the stress ensembles 
are probably interdependent [g. However, it seems that 
expectation values of quantities that depend on the ge¬ 
ometrical state only and not on the stress state are well 
described by the volume ensemble [H, possibly because 
the force-moment tensor can be treated as approximately 
constant in the systems considered here. Therefore we fo¬ 
cus on the volume ensemble in most of the present paper. 

Some doubts about the temperature-like interpretation 
of the compactivity have been raised on the basis of equi¬ 
libration experiments performed by Puckett and Daniels 
m- They found that in a two-component system, the 
compactivity did not equilibrate, whereas the angoricity 
did. There are different ways to explain this observation, 
some of which do not necessitate to give up the interpre¬ 
tation of the compactivity as a temperature-like variable 

m- 

A key assumption of standard statistical mechanics is 
the equivalence of time and ensemble averages, but me¬ 
chanically stable granular configurations are static states 
without any intrinsic time evolution. On the other hand, 
by applying the same external excitation to the granu¬ 
lar material again and again (i.e. tapping |S] or shearing 
[13] )■ this external excitation may take over the role of 
thermal agitation and the concept of a time average be¬ 
comes meaningful for granular statistics as well. Some 
tests of the ergodicity of granular systems are available 
in the literature. With systems of frictional discs excited 
with flow pulses, equivalence between time and ensemble 
averages was found M- In the case of a vertical tapping 
protocol, dependency on the protocol parameters was 
noted. Especially for small tapping amplitudes, noner- 
godicity was observed in numerical simulations |15j . This 
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may be related to the occurrence of irreversible branches 
in tapped granular systems [3 [50] . 

Several methods were proposed to determine the com- 
pactivity from experimental or simulation data [5J ll6H18j 
and applied to different kinds of granulates. In this work, 
we employ two-dimensional discrete-element (DEM) sim¬ 
ulations using polygonal particles to apply two differ¬ 
ent excitation protocols to otherwise identical granulates. 
This allows us to determine whether the calculated com- 
pactivity is independent of the specific excitation pro¬ 
tocol (which would support the Edwards theory) or not 
(which would oppose it). Recently, some authors have 
addressed the issue whether or not it is necessary to in¬ 
troduce protocol specihc extensions of Edwards-like ap¬ 
proaches mun]. We briefly discuss a possible approach 
in section IIB At least for the two protocols considered 
in this paper, we hnd protocol independence. Further¬ 
more, we test whether ideal-gas-like prediction analogies 
of the Edwards theory are consistent with the simulation 
data, and find good qualitative agreement, which even 
becomes quantitative if certain parameters are chosen by 
fitting. 

The paper is organized as follows. In sectionjnj we give 
an overview of some aspects of Edwards’ theory, relevant 
to this work. In section III we introduce the simulation 
technique. Section |IV| gives details of the applied exci¬ 
tation protocols and in section jVj the simulation results 


will be presented, evaluated and discussed. In Sec. VI 


some limitations of the volume ensemble are illustrated. 
Finally, section m gives a summary of our findings. 


II. THEORETICAL FRAMEWORK 


A. The microcanonical and the canonical volume 
ensemble 

The main assumption of the Edwards theory is the fol¬ 
lowing mm- If a granular ensemble is generated by a 
reproducible preparation protocol, all resulting mechani¬ 
cally stable configurations of the granular system which 
occupy the same volume will occur with the same prob¬ 
ability on repetition of the protocol. This means that in 
Edwards’ granular statistics the volume plays the same 
role as the energy in ordinary statistics. Consequently, 
the entropy S of the microcanonical ensemble having a 
certain volume V is given as the logarithm of the number 
n of stable states (in the permitted phase space) occu¬ 
pying this volume. 

Let a microstate of the granular ensemble, comprised 
of N particles, be described by a set of variables q. The 
entropy of this state is given by (see e. g. mmmy- 

S{V,N)=lnn{V,N), (I) 

n{V,N)= [ dq6{V-W{q)). (2) 

The integral runs over all stable conhgurations {< 7 } and 
the function W{q) gives the volume of the state q. Note 


that in this context, W{q) is the granular analogue of 
the Hamiltonian and V is the analogue of the internal 
energy. 

By analogy with standard statistical mechanics, an in¬ 
tensive, temperature-like variable x can be defined, usu¬ 
ally called compactivity: 


x = r' 


dV 

^ ■ 


( 3 ) 


In this paper, we will, for simplicity, mostly use the “ther¬ 
modynamic beta”, defined as the inverse compactivity 
P = 1 /x, instead of the compactivity itself. 

Note that in general, a constant A analogous to the 
Boltzmann constant may be introduced in Q and (|^. 
Here, we use A = I which means that we measure com¬ 
pactivity in units of volume. 

As in ordinary statistics/thermodynamics (see e.g. 
m) one may switch from the microcanonical ensemble 
to the canonical ensemble by a Legendre transformation. 
In the corresponding canonical ensemble, the probability 
of a microstate q occupying the volume W (q) is given by 
a Boltzmann-like distribution: 

P{q) = (4) 


with the canonical partition function 


Z = [ dq . (5) 

Once again, the integral runs over all possible mechani¬ 
cally stable states. Note that this is equivalent to a nota¬ 
tion often used in the literature, where the integral runs 
over all states and contains an additional factor 0{q), 
which takes the value zero for forbidden states and one 
for allowed ones, thus selecting the permitted states. 

Note that a tapping protocol does not necessarily lead 
to canonical sampling of the system. A trivial exam¬ 
ple would be tapping with so small an amplitude that 
the system is trapped in the initial blocked state. A 
non-trivial example is the irreversible branch for vertical 
tapping observed by Nowak et. al. [ 8 |. Also, the com¬ 
parison between the analytically solvable Bowles-Ashwin 
model system [5] and simulations where such a system is 
vertically tapped [lOj exhibits deviations from the canon¬ 
ical ensemble prediction. Note that the Bowles-Ashwin 
model assumes a highly confined geometry with much 
lower complexity than realistic granular systems. Ver¬ 
tical tapping applied to this special system with strong 
confinement may be unsuited for phase space exploration 
according to a flat probability measure. Or else the flat 
Edwards measure does not hold in general. Even if that 
were the case, a meaningful definition of compactivity 
might still be possible as will be described now. 


B. Generalised (protocol dependent) ensembles 

In case it turned out that Edwards’ assumption of a flat 
probability measure is not satisfied for the microcanical 
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ensemble, it is possible to modify the ensemble with (in 
general, state and protocol dependent) weighting factors 
w{q) [21 nni 1 ^ . The microcanical state density would 
be modified as follows 


nG{V,N)= f dqw{q)d{V-W{q)). (6) 
Providing we have 

^GiVi + V2,N) = naiVuNi) nGiV2, N 2 ) , (7) 


which is a much weaker assumption than a flat probabil¬ 
ity measure [3] , it is still possible to define a compactivity. 
The canonical formulation for the volume ensemble then 
reads 

P{q)=w{q)^e-^^^'iK ( 8 ) 

Zg = [ dq . (9) 

•7{9} 


A similar approach is feasible for the force-moment and 
the combined ensembles. Generalised canonical ensem¬ 
bles of this kind are used in statistical genetics [23l |24] . 

The qualitative impact of these modifications on the 
resulting thermodynamics is small. Especially equations 
0 and ( |2^ remain unaffected, as long as the weight¬ 
ing factors are the same for all data samples. If only one 
protocol is used, it is impossible to detect protocol de¬ 
pendencies in the weighting factors. But a comparison 
of different protocols, applied to the same system, may 
help to decide whether these factors, if their introduc¬ 
tion should turn out necessary, are protocol independent 
or not. 


C. Mean volume and volume fluctuations 


The volume W{q) of a certain state can be written as 
a sum of the mean volume V = (W) and the deviation 
from the mean: W = V + AP. Under the assumption 
that the volume fluctuations are small compared to the 
mean volume, we can write: 


U -b AU V U2 


(14) 


The mean volume fraction is therefore: 

^ =('/') = y , (15) 

, together with (AU^) = ay, we 

obtain 


and from eqs. ( 14|15 


(^1= ■ (16) 

Substituting Oy from eq. ( |l^ and reexpressing the dif¬ 
ferential of {V) according to d{V) = —Vg/{V)‘^d^, we 
find a relation between the mean volume fraction and its 
fluctuations: 
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<i> 


(j)^ d(j) 

VgTp- 


(17) 


Measuring the volume fraction fluctuations as a func¬ 
tion of the mean volume fraction and integrating equa¬ 
tion 0 or ( [I^ , respectively, is a way to calculate the 
compactivity up to an unknown constant. It is used fre¬ 
quently, e.g., in [SlITlld^. Note that determining the 
granular temperature via eq. (12) or eq. is just a rule 


of calculation, provided by Edwards’ theory but no proof 
of the theory. On the other hand, if after determining /3 
in some different way the relationship ( |I7| ) linking it with 
the volume fraction were not satisfied, a contradiction to 
Edwards’ theory would have been demonstrated. 


The calculation of the mean volume and its fluctua¬ 
tions is straightforward. The first derivative of the loga¬ 
rithm of ([^ (or (|^) with respect to (3 essentially is the 
mean volume 

^ ^ , ( 10 ) 

Op ^ J{q} 

whereas its second derivative is the variance of the vol¬ 
ume distribution, i.e., a measure for the strength of fluc¬ 
tuations. 

al. = (V‘)-{Vf = -^lnZ, ( 11 ) 

4 = -|(r). ( 12 ) 


Instead of the volume itself, we use the the volume frac¬ 
tion <f>, defined as the sum of the grain volumes Vg divided 
by the volume W occupied by the granulate: 


Hq) 


W{q) ■ 


D. Overlapping histogram method 

Another way to determine the inverse compactivity is 
the overlapping histogram method proposed in 2003 by 
Dean and Lefevre [25] ■ This method may also be used 
as a test whether or not a distribution of blocked states 
is Boltzmann-like distributed. In the original paper, the 
method was applied to the energy of the Sherrington- 
Kirkpatrick model for spin glasses, driven by a tapping¬ 
like mechanism, which has some similarities to granular 
dynamics. The method was then frequently used to de¬ 
termine the granular temperature mi in mi HH]- Under 
the assumption that Edwards’ theory holds, the proba¬ 
bility to measure a certain volume V in a granulate with 
an inverse compactivity /3q and a fixed number of grains 
N is 

P(V, ft,«) = 4, S(V- wm 

ZiPo,N) 


(13) 


(18) 
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where 53 (V, N) is the number of blocked states with vol¬ 
ume y, i.e S)(y,-/V) = fl{V,N). This equation holds 
even for the generalised non-flat probability measure ap¬ 
proach, when eq. §) is used for ^{V,N). 

We define the quantity Q as the logarithm of the ratio 
of probability densities for the volume V to arise, on the 
one hand at an inverse compactivity f3i and on the other 
hand, in a reference system held at a different inverse 
compactivity /3q: 


Q 

Q 


, P{V,PuN) 


= iPo - Pi)V + In 


:=Ai 


ZWo,N) 

ZWi,N) 


: = Bi, 


(19) 

( 20 ) 


Therefore, if we measure the probability density 
P{V,P,N) for different compactivities and calculate 
Q{V), the resulting function must be linear if Edwards’ 
theory holds. Evaluating the slope Aio allows us to deter¬ 
mine the inverse compactivity up to an additive constant. 
Furthermore, the intercept Bio is nothing else than the 
logarithm of the partition function up to an additive con¬ 
stant at the corresponding inverse compactivity. This 
permits testing the validity of Edwards’ assumption in 
the sense, that if Q{V) were not a straight line, neither 
Eq.Q nor Eq. would describe the probability density 
of the system correctly. 

However, one has to be somewhat careful with the in¬ 
terpretation of the results, as has been pointed out in 
m- Under certain circumstances, very similar results 
as the ones expected from Edwards’ theory can occur, 
if the distributions of the samples are just Gaussian. In 
appendix [5] we discuss this situation. 


E. Ideal quadron solution 

There are very few real ab initio predictions from Ed¬ 
wards’ theory in the literature [6l |22l [29l [30], due to 
some general difficulties. In order to calculate analyti¬ 
cal expressions for the partition function, knowledge of 
an explicit expression for the granular Hamiltonian W (q) 
is necessary. Of course, if all positions, orientations and 
shapes of the grains are known, the occupied volume is 
a function of these quantities, but in practice it is not 
easy to write down an explicit equation and even if this 
can be done, the integration over the permitted blocked 
states is very difficult because the permitted states are 
unknown in general. 

Therefore, attempts to calculate the partition function 
were based on standard volume tessellations, such as the 
Voronoi and Delaunay tessellations [5DH55] . A possible 
alternative is an arch-based approach m which a priori 
takes only stable configurations into account. Blumenfeld 
and Edwards proposed a physically motivated tessella¬ 
tion based on the quadron construction [SmilMlIM]- In 
principle, quadrons are used as quasi-particles describing 


the structure of the granulate at any arbitrary position 
within the system in a distinct way. It was mentioned 
in the literature that the ideal quadron tesselation fails 
in the presence of non-convex voids in the granulate [35j . 
However, in a system of monodisperse spheres such non- 
convex voids vanish by neglecting rattlers [36] . Even if 
some non-convex voids remain, they could be tesselated 
by convex polygons which repairs the quadron tessella¬ 
tion, at the price that the number of quadrons increases. 
As long as non-convex voids are the exception rather than 
the rule, they will not produce significant changes to the 
calculations. 

Under the (very rough) assumption that quadrons oc¬ 
cupy volumes between Vo~ A and Vq-I-A at constant den¬ 
sity of states and that there are no interactions between 
the quadrons, the partition function can be calculated 
explicitly for two dimensional systems (see [29]): 




( 21 ) 


where N is the number of particles and z is the mean 
coordination number in the granular system. This ap¬ 
proximation is called the ideal quadron approximation 
by analogy with the description of ideal gases in ordi¬ 
nary statistics. Note that the partition function (211 is a 
special version of a more general ideal-gas-like approach. 
If one assumes that the volume is tesselated by a number 
N of statistically independent, non-interacting elemen¬ 
tary cells and their volume is restricted to an interval be¬ 
tween a minimal volume Vq — A and a maximal volume 
Vo -|- A, without any additional assumption on the na¬ 
ture of these elementary cells one ends up with Eq. (211, 
on replacing Nz^N. Contrary to the very general 
ideal-gas-like approach, there are some possibilities to go 
beyond the interaction-free situation in the quadron ap¬ 
proach, which will be considered in future work. 

Using (10), the ideal quadron prediction for the mean 
volume is obtained: 


(U) = Nz[Vq + -- Acoth(,dA) 


( 22 ) 


For the current work, it is helpful to rewrite (22) in terms 


of the volume fraction. Dividing (22) by the total grain 


volume Vg results in 


- 1 NzVq Nz NzA , 

- -f —-— coth(;3A) . (23) 




VgP 


Va 


In order to specify the free parameters Vq and A, we as¬ 
sume that the limits of can be identified with the 
volume fractions of random loose packing (rip) and ran¬ 
dom close packing (rep) in the following way: 



lim (j) ^ 
/ 3^0 


lim (j) ^ 
^—¥00 


NzVo 

Nz{Vo - A) 


(24) 

(25) 
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Now we can rewrite (231 as 


A-l 


+ § - N' 


(26) 


where Below, eq. p6|) will be useful 


as a fitting function for our simulation data. 

We remark that the ideal quadron solution is a very 
rough estimation, because the main difficulty in calcula¬ 
tion of the partition function ([^ is filtering out stable 
configurations. In the ideal quadron solution, a minimal 
filtering is performed in the sense that there is a minimal 
and a maximal volume per quadron so that configura¬ 
tions that are either too loose or too dense are filtered 
out. However, the number of states between these lim¬ 
its may be overestimated. The arch-based approach m 
may have the capability to overcome this issue. However, 
in the current form no analytical equation such as (26) is 


available. Only numerical solutions under very simplify¬ 
ing assumptions are on hand, which makes a comparison 
of our data with this model difficult. 

Note that there is no commonly accepted definition 
of the states rip and rep. In this paper we use these 


terms in the sense of eq. (241 and (25). An interesting 


fact is that in the framework of the ideal quadron model 
for /3 —> 0, which is the limit of very high compactivi- 
ties, the mean volume per quadron is Vq , not Vq + A. 
States which would have a mean volume per quadron 
with V > Vo correspond to a population inversion, i.e., 
states with negative granular temperature. It may be 
speculated that these states correlate with so-called ran¬ 
dom very loose packings states [SS] , which are states only 
achievable with very special protocols. A verification of 
this idea is beyond the scope of the present article. 


III. SIMULATION METHOD 

In our simulation, we extend an existing DEM Code, 
originally developed for the investigation of the me¬ 
chanical properties of granular piles consisting of two- 
dimensional polygonal particles USED]. The dynamics 
of the Ah particle’s position and orientation (j)i are 
described by Newton’s and Euler’s equations of motion: 



FIG. 1. (Color online) Sketch of two polygonal particles in 
contact. The normal and tangential direction of a collision is 
determined by the contact line (green), the contact point is 
defined as the middle of the contact line. 


For fast determination of potential particle contacts, 
the particles are surrounded by bounding boxes and 
we employ an incremental sort-and-update algorithm 
ESES] to identify overlapping bounding boxes efficiently. 
Whenever the bounding boxes of two particles overlap, 
we use a closest-feature algorithm ESEDj from virtual 
reality and robotics applications jUj to calculate the 
polygon distance. In the worst-case scenario, the com¬ 
putational complexity is O(nlogn) where n is the num¬ 
ber of polygon features (corners and edges). However, 
the typical behavior, whenever a good guess from the 
last timestep is available, is the calculation of the poly¬ 
gon distance in constant time 0(1), independently of the 
polygon edge number. 

Figure shows a sketch of two particles in contact. 
The contact point Sij between two particles is defined as 
the midpoint of the line between ci and C 2 . For the force 
calculation, we define some quantities first: the charac¬ 
teristic length 


I = 


ri + Tj 


(29) 


(note that I is not the length of the contact line), the 
reduced mass 


m± = -, 

rUi + rtij 


(30) 


^ , (27) 

Jii'i = X! • ( 28 ) 

Herein, Fij is the contact force between particle i and 
particle j, is the corresponding torque acting on the 
particle (referred to its center of mass), and F^ is the 
external force acting on the particle (e.g., gravity). The 
mass of particle i is denoted by rrii and its moment of 
inertia by Ji ETHID] . We assume external torques to be 
absent. 


and the reduced tangential mass, including the moments 
of inertia 


— + — + — -I- — 

rrii rrij Ji Jj 

The relative tangential velocity at the contact point Sij 
is 

U|| = (v^ - Vj + (vi X (jJi) - {Vj X Ljj) ) • n||. (32) 

where v^, Vj are the velocities of the particles and ujj 
are their angular velocities. Furthermore, we define the 













effective penetration length as 


IV. EXCITATION PROTOCOLS AND 
PARAMETERS 
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/leff = y . (33) 

Note that I does not change significantly during a col¬ 
lision, so hgff is essentially proportional to the overlap 
area A. The normal component of the contact force is 
determined by the equation 


F± = max 


Ehes - 'y\/Em ±heff 


0 


(34) 


Here E is the two-dimensional Young’s modulus and 7 is 
the dissipation strength. 


The maximum function in (34) ensures that the nor¬ 


mal force cannot become attractive. Physically, situa¬ 
tions where the dissipative term overcomes the repulsive 
force correpond to the case that the particles’ separa¬ 
tion velocity is faster than the relaxation velocity of the 
grain deformation, i.e., the contact between the particles 
is lost, before the overlap of the non-deformable model 
particles becomes zero [IS]. Consequently, the contact 
force vanishes from this moment on. 

We used the Cundall-Strack model [IS] for modelling 
tangential forces. At the beginning of the collision, the 
Cundall-Strack force E* is zero. If this force is given at 
the previous timestep, corresponding to time tc — At, the 
Cundall-Strack force at the current timestep at time tc 
is determined by 


F\\ (tc) = 


F|| (tc - At) -h n||(tc - At)At-A ; nFj_{tc) 


(35) 


One may visualize this model as a spring between the 
contact points being established when a new contact ap¬ 
pears and the length of the spring being limited by the 
value of Coloumb friction (Coulomb condition). If this 
value is reached, the points of attachment of the spring 
start to slide, so that the spring length remains constant. 
In this way, the Cundal-Strack model mimics sticking 
and sliding friction. In order to avoid unrealistic oscil¬ 
lations, a damping term proportional to the velocity is 
added to the Cundall-Strack force, which also satisfies 
the Coloumb conditions. The complete trangential force 
then is 


F||(tc) = nrin 


A'lj' (tc) + ^11 



t^F±{tc) 


(36) 


With (34) and 


all the contact forces and the re¬ 


sulting torques between the particles and between parti¬ 
cles and walls (which are treated as particles with fixed 
position and orientation) are defined. The equations of 
motion (27) and ([2^ are then solved numerically using 


a sixth-order Gear predictor-corrector method un¬ 


inspired by former experimental studies |I6H18] . we 
implemented two protocols for exciting the granular mat¬ 
ter periodically. In both protocols, an excitation period 
in which the phase space is explored alternates with a 
relaxation period in which the grains come to rest com¬ 
pletely. Both protocols are characterized by a control 
parameter, which we call tapping parameter. 

In all simulations, we used a bidisperse mixture of 1184 
regular decagons, composed of 544 particles with radius 
ri = 9 mm and 640 particles with radius r 2 = 6.36 mm. 
Young’s modulus was set to E = 1000 N/m and the sim¬ 
ulation time step was chosen as At = 5 • 10“^ s. The 
Coulomb friction coefhcent was taken to be /r = 0.6 and 
for the normal friction coefhcent 7 = 0.5 was used. The 
particles mass density was p = 0.001 g/mm^ and was the 
same for big and small particles. 

With both protocols, the width of the simulation area 
was 280 mm. In the case of the negative g protocol, there 
is no lid, for the rotating g protocol the height of the box 
is 900 mm. Walls are realised as hxed rectangular parti¬ 
cles with the same Young’s modulus and friction coefh- 
cients as the mobile particles. 

Tejada et. al. [48] pointed out that the size of the time 
step in DEM simulations may influence the width and the 
shape (but not the mean) of the determined probability 
distributions, even if the time step appears sufhciently 
small using common criteria. In order to make sure that 
this effect does not influence our results, we repeated 
some of the simulations with a time step of At = 10“® s. 
We found that the mean volume fraction, the volume 
fraction variance and also the shape of the volume frac¬ 
tion distribution did not change on reduction of the time 
step. 

The first protocol, called “pulses of negative gravita¬ 
tion” or “negative g” protocol, is illustrated in figure]^. 
Most of the time, the granulate is at rest in a container 
under normal gravitation, but for short time intervals, 
the direction of gravitational acceleration is reversed. 
The time dependence of the gravitational acceleration 
g{t) is taken to be 


g{t) = 9.81 m/s^ey • 

Here t = t mod T is the simulation time t modulus the 
period T of the protocol and ti is the duration of a pulse. 
We took the pulse length to be ti = O.I s and the period 
T = 3s. The relaxation period was chosen so that for 
the biggest excitation amplitudes the resulting volume 
fraction variations due to remaining kinetic energy in the 
system are smaller than 10 “^ which is 100 times smaller 
than the typical magnitude of the observed static volume 
fraction fluctuations. The pulse amplitude gp was used 
as a tapping parameter. Due to Einstein’s equivalence 
principle, this protocol is equivalent, per period, to a 


t < ti 
otherwise 


(37) 
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(b) 





27r/a; 


t 


FIG. 2. (Color online) a) The “negative g” protocol: For a 
short time interval ti the gravity is turned upward with a 
prescribed magnitude. Afterwards, the granular system has 
time to come to rest completely, before the next pulse starts, 
b) The “rotating g” protocol: The gravitational acceleration 
vector performs a complete rotation. Afterwards, the granular 
system comes to rest completely before the next round starts. 


downward acceleration {gp + l)g for an interval ti and a 
subsequent stoppage in the gravitational field for a time 
span T—ti, taken long enough for the granulate to come 
to rest. Here and in the following when we write g, we 
mean g = 9.81 m/s^. 

Figure illustrates the second protocol. The gran¬ 
ulate is at rest in a closed box. Then the gravitational 
acceleration performs a complete rotation followed by a 
relaxation period. 


g{t) = 9.81m/s^ 


sin(wt)ea; — cos{ujt)ey i < ^ 
—By otherwise 


(38) 

As in the negative g case, i = t mod T is the simula¬ 
tion time t modulus the period T. The angular fre¬ 
quency uj is the control parameter of the protocol. Of 
course T > 2tt/uj must be fulfilled. In our simulations, 
we used T = 27r/a; -I- 3 s, with ut > 4s“^. For both 
protocols, we tested whether or not segregation of the 
particles occurs by measuring the cumulative number of 
small and big grains as function of the column height 
n{h) (i.e. the number of small or big particles with z co¬ 
ordinates smaller than h). These curves are straight lines 


and do not change during the simulations, except for fluc¬ 
tuations. Therefore we conclude that in our simulations 
segregation did not take place. 


V. DETERMINATION OF COMPACTIVITY, 
FLUCTUATIONS AND THE PARTITION 
FUNCTION 

For both protocols and each choice of the tapping pa¬ 
rameter the simulations ran for 2500 — 3000 excitation 
and relaxation periods. Note that this relatively high 
number of taps is necessary to draw serious conclusions 
for the systems used here. While the confidence interval 
of the mean volume fraction became sufficiently small 
after some hundred taps, this was not the case for the 
estimated variance. In test simulations, where only a 
few hundred taps were considered, the uncertainty of the 
estimated variance [SSJ pp. 771-772] was as big as the 
domain of the measured variance itself. 

Immediately before the relaxation period ended, we 
determined the volume fraction of the granulate by mea¬ 
suring the fraction of solid particles in a test volume. The 
test volume shown in figure]^ is a square with edge length 
of 400 mm and was chosen in such a way that some layers 
of particles were between the borders of the test region 
and the walls. 



FIG. 3. (Golor online) A typical situation when the granulate 
is at rest. The region shaded grey is the test volume, used for 
volume fraction determination. 

It was pointed out [IH] that test volumes must be big 
enough to avoid size depended effects. To make sure that 
our test volume is sufficiently large, we divided it into two 
neighbouring columns. The relative difference of mean 
volume fraction (after 2500 taps) between the columns 
and the entire test volume was always smaller than 0.01 % 
and the deviation of the volume fraction fluctuation ratio 
from y/2 was always lower than 1 %. 

Figure shows some typical time series for both pro¬ 
tocols. The mean volume fraction reaches a steady state 
value very quickly (after approximately < 10 taps) and 
then only fluctuations around its mean value occur. 
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FIG. 4. (Color online) Exemplary time series for both pro¬ 
tocols. The points are the volume fractions calculated from 
simulation data and the solid lines are the corresponding mean 
values. The tapping parameter Is for the “negative g” pro¬ 
tocol (a): Qpg — 6m/s^ (red points) and gpQ = 28m/s^ (blue 
points) and for the “rotating y” protocol lj = 0.75 • 27r s“^ 
(red points) and tj = 2 • 27r s“'^ (blue points). 


In Fig. the mean value of the volume fraction 0 
and the standard deviation of the volume fraction dis¬ 
tribution, characterizing the volume fraction fluctuation 
strength, are shown. Initially, the volume fraction de¬ 
creases with increasing pulse strength until it reaches a 
local minimum around (?p = 20. Afterwards, the volume 
fraction increases with the pulse strength. Similar be¬ 
haviour for tapped granulates at high tapping strengths 
was described in isn Eg. A possible explanation is 
the interplay between two competing effects. First, the 
stronger the pulse of negative gravitation, the more the 
granulate is whirled around, i.e. the looser the resulting 
packing should get. On the other hand, the stronger the 
pulse, the higher the granulate flies, therefore the higher 
its impact velocity when it reaches the bottom, result¬ 
ing in stronger compaction during the relaxation period. 
The first effect dominates for relatively small values of gp, 
the second effect becomes more important for stronger 
pulses. At the same value of gp, where the volume frac- 


FIG. 5. Simulation results for the “negative g” protocol: (a) 
the mean volume fraction as function of the pulse strength in 
units of earth acceleration is depicted, (b) shows the standard 
deviation of the volume fraction fluctuations. The error bars 
correspond to a confidence interval of 95 percent. 


tion is minimal, the volume fraction fluctuations have a 
local maximum (see Fig. Ef)- Similar results were found 
in work by Pugnaloni et al. |52j . 

The corresponding results for the “rotating g” proto¬ 
col are shown in figure With increasing frequency the 
mean volume fraction shows a crossover from 0 ~ 0.814 
to 0 ~ 0.821. When the frequency is very small, a com¬ 
plete rotation takes a long time and the granulate be¬ 
haves in a quasistatic way. It is clear that in this case a 
further decrease of the frequency will not have a signif¬ 
icant influence on the resulting volume fraction. There¬ 
fore, for small frequencys the 4>{ui) curve should approach 
a plateau. When the frequency increases, the rotation 
gets faster, the granulate is more strongly whirled around 
and becomes looser. When the frequency increases fur¬ 
ther, the time per rotation gets smaller which compen¬ 
sates the increase of swirling due to faster rotation. For 
very high frequency, one expects that the system reaches 
an irreversible regime, because it does not have time to 
respond to the rotation pulse and will largely stay in the 
initial configuration - but this regime is not reached in 
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FIG. 6. Simulation results for the “rotating g” protocol: Part 
(a) shows the mean volume fraction as function of rotation 
frequency / = uj /{ 2 - k ) and (b) shows the standard derivation 
of the volume fraction fluctuations. The errorbars correspond 
to a confidence interval of 95 percent. 


our simulations. 

In order to determine the granular temperature of 
the samples with the overlapping histogram method, the 
probability density distribution of the volume fraction 
must be estimated. Notwithstanding the name of the 
approach, we used the more sophisticated kernel density 
estimation method [531154j instead of histograms in order 
to obtain an approximation for the probability density. 
A normal kernel was employed and the bandwidth was 
chosen according to Silverman’s rule of thumb |55| . In 
appendixj^ a short description of the approach is offered. 

Some of the determined probability density distribu¬ 
tions are shown in Fig. for the “negative g” protocol 
and in Fig. for the “rotating g” protocol. In order 
to test if the distributions are Gaussian, we used a chi- 
square test [35] with a significance level of 5%. The null 
hypothesis that the data is normally distributed was re¬ 
jected for all samples m- Although Gaussians may still 
be good approximations to the central region of our dis¬ 
tributions, we take this as evidence for a statistical me¬ 


chanical origin of these distributions rather than their 
emergence from some unknown additive process not re¬ 
lated to phase space exploration. Therefore, we believe 
that our data do not give spurious results due to the 
pitfall described in appendix [A| 

We also show in appendix^ that while the presence 
of Gaussian distributions may lead to false positive re¬ 
sults regarding the validity of Edwards’ theory this 
is by no means true for all Gaussian distributions: in the 
limit of large numbers, the canonical distribution corre¬ 
sponding to a fixed compactivity must become Gaussian 
as well, and this distribution obviously satisfies Edwards’ 
theory by definition. We take the fact that the tails of 
our distributions are non-Gaussian together with the ver¬ 
ification of (201 as an indication that the theory works. 


This interpretation is supported by other studies of very 
similar systems where the volume-per-particle distribu¬ 
tion (which is not a sum and therefore not subject to 
the central limit theorem) was found not to be Gaussian 
[IHIIIT] even in the central region. 


In terms of volume fraction, equation (201 reads: 


Q = In 


m/^i) 

Pi^,M 


fo 1 ^(/3i) 


(39) 


We note that (391 is interpreted in terms of a canoni¬ 


cal ensemble, which implies the number of particles to 
be a constant. This means that the total volume of the 
grains Vg in the test volume must be a constant jSS]. In 
principle, one would have to adapt the test volume size 
Vt used in determining the volume fraction distribution 
to the measured mean volume fraction (j) in such a way 
that Vg = Vt y remains constant. The size of the test 
volume would then be a function of the mean volume 
fraction. This could become important, if very small test 
volumes were used. The standard deviation of the vol¬ 
ume fraction distribution should decrease with increas¬ 
ing test volume as cr oc which follows from ( |T7| , 

while (j) does not depend on the system size. By using 
a constant test volume VV, the magnitude of differences 
in Vg due to different volume fractions is bounded by 
AEg = Er(l/^('rcp — l/'/’rip) — O.OSEr- If the test volume 
is large compared to the size of the particles, the rela¬ 
tive error which is made by using a constant test volume, 
assuming that the cumulative grain volume is constant 
(in spite of the different volume fractions), corresponds 
to approximately (2 — 3)%. This is smaller than the con¬ 
fidence interval of tr^ due to the finite sample size and 
therefore negligible. 

By choosing a reference probability distribution which 
is used as the denominator in (39) we are able to de¬ 


termine the inverse compactivity and the logarithm of 
the partition function up to additive constants, which 
are the unknown inverse compactivity of the reference 
distribution and the unknown logarithm of the partition 
function thereof, respectively. In order to avoid errors 
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FIG. 7. (Color online) (a) Probability density function of 
the volume fraction for the “negative g” protocol, estimated 
with the kernel density method, for different values of the 
pulse strength, (b) The quantity Q from equation (201 for 
some values of the tapping parameter, where the distribution 
with gp = 6 was chosen as denominator in (19l. (c) The 

compactivity calculated with the help of the overlapping his¬ 
togram method using the distribution with gp = 6 (green, 
squares) and Qp — 10 (red, circles) as denominator in (|19[). 


due to insufficient data in the tails of the distribution, 
we evaluated the slope of Q only in regions where the 
value of each distribution involved in the calculation of 
( [M| ) is bigger than 5% of its maximum. 

Figures and[^ show the quantity Q for some sam¬ 
ples as function of the inverse volume fraction, where the 
distribution with gp = 6 oi the “negative g” protocol was 
chosen as the reference distribution, because it has suf¬ 
ficient overlap with all the other distributions. We tried 
several other distributions as references, too. So long as 
the distribution overlap was big enough, we always found 
a linear relation between Q and the inverse volume frac¬ 
tion 4’~^. This may be interpreted as suggestive of the 
validity of Edwards’ assumptions. Note that the devia¬ 
tions from the straight line on the left and right ends of 
the curves are due to insufficient number of sampling data 
points in that region. Therefore the estimated probabil¬ 


FIG. 8. (Color online) (a) Probability density function of the 
volume fraction for samples of the “rotating g” protocol, es¬ 
timated with the kernel density method, for different values 
of the rotation frequency, (b) The quantity Q from equa¬ 
tion ( |20[ ) for some values of the tapping parameter using the 
distribution with gp = 6 of the “negative g” protocol as de¬ 
nominator in (19l. (c) The compactivity calculated with the 
help of the overlapping histogram method using the distribu¬ 
tion from the “rotating g“ protocol with / = 0.5Hz (blue, 
circles), / = 1.25Hz (green, squares) and the distribution 
from the “negative g” protocol with gp = 6 (black, diamonds) 
as denominator. 


ity density function behaves like the tails of the sampling 
kernel, which is reflected also in Q. This behaviour is 
a mathematical artefact and not a systematic deviation 
from a straight line. 

From (39), it follows immediately that the values de¬ 
termined for /3 using different samples as reference distri¬ 
bution may only differ by an additive constant. This pre¬ 
diction holds for the “negative g” protocol as is demon¬ 
strated in Figure!^. There, the samples with gp = 6 and 
gp = 10 were usea as reference distributions. The same 
is true for the “rotating g” protocol (Fig. |^) when we 
used samples corresponding to different values of the ro¬ 
tating frequency. It even applies if we use samples from 
the “negative g” protocol as reference distribution, in 
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FIG. 9. (Color online) The compactivity as function of the 
volume fraction density for the “negative g” protocol (square) 
and for the “rotating g” protocol(diamonds). The values were 
corrected with the additive constant Po determined from the 
ideal quadron fit. The solid line is a fit to the ideal quadron 
solution 


agreement with the fact that the predictions of Edwards’ 
theory should be protocol independent. 

From now on, we use the distribution from the “nega¬ 
tive g” protocol with gp = 6 as reference distribution for 
all the following calculations. Figure]^ shows the values 
determined for the inverse compactivity for both proto¬ 
cols as function of the volume fraction. All the data, 
whether obtained from the branch left of the minimum 
or the branch right of the minimum in the 4>{gp) curve 
of the “negative g” protocol (cf. Fig. or else from the 
“rotating g” protocol is fitted by the same function /3((()). 
This is a strong indicator for the applicability of Edwards’ 
theory to our samples. 

In order to ht our results to an analytical function, we 
used the ideal quadron solution ( [^ . The values of ran¬ 
dom loose packing (rip) and random close packing (rep) 
were assumed to be (pj-ip = 0.81 and (py-cp = 0.855. The 
choice of these values was pragmatic. To our knowledge, 
there are no studies about the exact values for random 
close and random loose packing for bidisperse decagons. 
Furthermore the values of rip and rep will depend on on 
the size ratio between the particles. Therefore, we used 
plausible values obtained for (bidiperse) disks [13 EH ED]- 
We checked that the influence of the values of random 
loose and random close packings on the values of the ob¬ 
tained parameters are smaller than 20%, as long as the 
values are in the plausible interval. 


Nz/Vg was used as htting parameter. Since we can 
determine the inverse compactivity only up to an addi¬ 
tive constant from the overlapping histogram method, we 
made the replacement /3 —>■ /3 -I- ,5o in (26), where $ is the 
value determined from the simulations and /3o is taken 
as an additional fit parameter. Note that the determi¬ 
nation of Po is possible as we assumed that the inverse 


compactivity is fixed at random loose and random close 
packing and that this is not specihe to the ideal quadron 
fit. We emphasize that the parameter Pq only shifts the 
whole curve shown in Fig. [^ upward or downward and is 
the same for both protocols. Since the same value of Pq 
was added to every data point obtained from the over¬ 
lapping histogram method for the comparison of the fit 
function with the simulation data in Fig. this value 
does not influence the conclusion that the compactivity 
is the same for both protocols. The fit describes the sim¬ 
ulation data very well, as is seen in Fig. However, 
we find as htting value Nz/Vg = 556.51m , which dif¬ 

fers by about a factor of 40 from the values used in the 
simulation. This might be understood by assuming that 
approximately 40 quadrons constitute a statistically in¬ 
dependent unit in the granular ensemble. This issue cer¬ 
tainly needs further study. We also tried to Ht pyip and 
(/)rcp and using the known particle number in (261 but this 
leads to an unphysical value for random close packing of 
about 4>rcp — 1-2. 

When we used only the data obtained by one of 
the protocols to determine the parameters of the ideal 
quadron solution, the solution Htted the data of the other 
protocol. The relative deviation of the obtained Htting 
parameters is smaller the 2 %. Because this results in 
curves that are indistinguishable to the eye, only the Ht 
using the whole dataset is presented in in Fig. Never¬ 
theless, this means that we can predict the P{(f>) curve, 
for the “rotating g” protocol using the data determined 
from the “negative g” protocol and vice versa. However, 
for all compactivities determined, the same reference dis¬ 
tribution was used so the data of the Ht employed for the 
one and the other protocol was not entirely independent. 

While the slope of Q allows us to determine the in¬ 
verse compactivity, the axis intercept B determines the 
logarithm of the partition function 


B = -\nZ{P)+hiZ{po). 


(40) 


If the assumptions leading to ( |39| ) are correct, the par¬ 
tition function of the ideal quadron solution (21) should 


describe the found intercept without additional Htting. 
The parameters A, Vq, Nz are directly related to the pa¬ 
rameters determined via Eqs. (24) and (25). As Fig. 10 


shows, the numerically determined values and the ideal 
quadron solution Ht very well, independently of the pro¬ 
tocol and of the branch in the “negative g” protocol. 

Using (26) in ( [T7| ), we get the relationship between the 
mean volume fraction and its fluctuations. The result 
together with simulation data obtained directly is shown 
in Hgure [m The data is in good agreement with the 
theory for both protocols. 

We mention that in a former study which used verti¬ 
cally tapped monodisperse regular polygons m , a max¬ 
imum in the p — a curve was reported which coincided 
with an inflexion point in the impulse strength - volume 
fraction curve. In our simulation we do not see this ef¬ 
fect, also in experimental work about bidisperse two di¬ 
mensional systems such a maximum was not observed 
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FIG. 10. (Color online) The logarithm of the parti¬ 
tion function, determined from the overlapping histogram 
method(symbols). The solid line is the ideal quadron solu¬ 
tion. 


moment and volume ensembles (see, e.g. IS])> deviations 
may occur. 

VI. LIMITATIONS OF THE VOLUME 
ENSEMBLE 

Whereas the volume ensemble appears to succeed in 
describing the geometrical and structural degrees of free¬ 
dom of a granular aggregate, this is not the case for the 
stress state of the latter. If the volume ensemble entailed 
a complete description of a granular state and its prob¬ 
ability distribution, the mean stress of the system would 
have to be a unique function of the inverse compactivity 
and therefore also a unique function of the mean vol¬ 
ume fraction. We computed the mean extensive stress 
(or force-moment tensor), defined as 

( 41 ) 

p 



0 


as a function of the volume fraction. Note that the vol¬ 
ume density of this tensor is the stress itself. Here i and 
j label Cartesian coordinates. The sum runs over all par¬ 
ticles, where cr? is the mean stress in particle p and 
is the volume associated with the particle. The result is 
shown in Fig. The stress tensor is obviously not a 
unique function of the volume fraction. 


a 






0 0 


FIG. 11. (Color online) The volume fraction fluctuation (stan¬ 
dard deviation) as function of the mean volume fraction. The 
blue points correspond to the “rotating g“ protocol and the 
magenta points and lines corresponds to the negative g” pro¬ 
tocol. The solid line is the ideal quadron fit. 


[18]. It might be speculated that crystallization effects 
that occur frequently in two-dimensional monodisperse 
systems were responsible for the occurrence of the maxi¬ 
mum in Ref. m, but this question cannot be clarified in 
this study. 

However, we cannot exclude that there may be a small 
hysteresis for the branch pieces to the left of minimum 
and to the right of minimum, respectively (cf. Fig. [^. 
Clearly, if the volume ensemble completely described all 
structural degrees of freedom and the probability, two 
states with the same /3 and the same (/> would be iden¬ 
tical and therefore cr^ would also have to be the same. 
However, if the volume ensemble is only a good approxi¬ 
mation of the geometric aspects of interdependent force- 
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FIG. 12. (Color online) The components of the extensive 
stress tensor as function of the mean volume fraction for the 
“negative g” protocol (diamonds) and the “rotating g” pro¬ 
tocol (crosses). 

Neither the results from the different protocols nor the 
results from the half-branches left to the minimum and 
right to the minimum, respectively, of the “negative g” 
protocol fall on the same curve, which is in agreement 
with similar findings on tapped granular matter [52j . 
Two systems with almost identical particle positions and 
orientations can be in very different stress states which 
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is not captured by the volume ensemble. To describe the 
stress state and the volume state together, probably the 
combined volume-stress ensemble [B] must be taken into 
account. We remark that it is unclear so far wether or 
not the stress states observed here are “very different” 
or even “very similar” because we do not know the size 
of the accessible phase space for the extensive stress ten¬ 
sor. Maybe the deviations in our systems which have 
magnitude on the order of 0.01 Nm are so small that it 
is reasonable to assume in first approximation that the 
stress-moment tensor is almost constant which would be 
a possible explanation for the success of assuming a pure 
volume ensemble. 


VII. CONCLUSIONS 

We used two different protocols to excite a granular 
ensemble periodically. The inverse compactivity was de¬ 
termined as a function of the mean volume fraction and 
we found that the relation between the compactivity and 
the mean volume fraction is protocol independent. We 
determined an expression for the logarithm of the par¬ 
tition function and thus the thermodynamic potential 
which is the equivalent of the free energy of classical 
statistics. This was done by using the ideal quadron solu¬ 
tion derived by Blumenfeld and Edwards [5^ as a fitting 
function. Even though the ideal quadron solution makes 
very rough assumptions, the resulting description is in 
qualitative agreement with the findings from the simula¬ 
tions. If the particle number is used as a fit parameter 
instead of calculating the distribution with the true par¬ 
ticle number, the ideal quadron solution is also able to 
describe the results quantitatively. 

We found that all our simulation results related to 
structural quantities are compatible with Edwards’ the¬ 
ory and that the Edwards theory describes the volume 
fluctuations very well, independently of the excitation 
protocol. The usage of the Edwards volume ensemble 
seems to be sufficient for the description of system prop¬ 
erties which are related to the geometrical arrangement 
of the grains, but as might be expected from previous 
findings, it is not sufficient to describe the stress states 
of the granular ensemble. 


Appendix A: Overlapping histograms for Gaussian 
distributions 


If two volume samples were Gaussian distributed with 
mean values Vi and V 2 and variances af and cr|, the 
corresponding Q-function defined in (|I^ would be [27]: 


Q^(V) = 


(V-Vi)^ 

2af 


ct" - V2? 
2 a| 



This is a quadratic function of U, but in some Y interval 
the curvature of Q® may be very small and the parabolic 


function (A1) would then practically be indistinguishable 


from a linear function. This happens in particular, if the 
variances G\ and cr2 are close to each other, i.e., for nearby 
compactivities. If we define a function as the slope 
of (A1) midway between the maximum values of the two 
normal distributions. 


^21 


we obtain EH; 




V={Vt+V2)/2 


(A2) 




Identifying formally 


= /3(^2) - 


(A3) 


(A4) 


we find from (A3): 


=21^ + ^ 
P{V2) - /3(IA) I fT? ^ a? 


-1 


(AS) 


If we assume that the variance is a unique function of the 
mean volume, this equation is an approximation of (12) 


which is the better the smaller the difference IV2 — Vi\7 
We note that the formal identification (A4) is strictly 


speaking inherently contradictory, which is easy to see if 
we calculate Q for a third sample with mean value V3 
and variance <73 and the same reference sample in the 
denominator. From ( |A4| ) A31 — A21 = A32 f ollows, but 
if one calculates the same quantity from (A3), the result 


does not agree. (However, the identification is possible, 
if the three variances are the same.) 

Analogously, we can calculate the intercept of the 
tangent which touches (A1) at U = (U1 -|- V2)/2: 






In — 


By calculating the limit 


R9 

lim -Al 


= -Vi - 


(9cri 


(A6) 


(A7) 


we find that the term is an approximation for 

(101, if we assume to be negligibly small and formally 


identify = ln(Z 2 ) — In(Zi). Again the approximation 
becomes better as the difference between the volume V2 
and the reference volume Vi gets smaller. 

Due to these similarities, it might occur that even if the 
relations (101, (12) and (20) are consistent, that within 


the limits of data precision it is not possible to decide 
whether the reason of this agreement is the correctness 
of Edwards’ theory or the fact that the data generated 
by a specific protocol happen to have a distribution that 
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is well approximated by a normal distribution. If the 
samples are generated using the same protocol, one may 
expect that there is a function but when different 

protocols are used, it would be surprising, if both proto¬ 
cols led to the same relation, unless a general principle, 
such as Edwards’ theory, were at work. 

On one hand, these similarities make it more difficult 
to verify Edwards’ theory, on the other hand, due to the 
central limit theorem, we should expect that the distribu¬ 


tion of Eq. (18) becomes more and more Gaussian with 
increasing system size. Therefore, it is not always true 
that the appearance of Gaussian distributions signifies in¬ 
applicability of the overlapping histogram method in the 
determination of the compactivity and of related quanti¬ 
ties. Let us briefly have a look at this. Rewriting Eq. (18) 


for general {3 and using the microcanonical result for the 
density of states, we have 

^ ^ Z{/3,N) Z{I3,N) ■ ^ ^ 

As in standard statistical mechanics, we can then argue 
that for large systems this distribution has a sharp peak 
at the mean value of the volume and we may expand the 
entropy about this average, neglecting terms that are of 
higher than quadratic order: 

^(E, N)-(3V^ S{V, N) + P{V) {V - V) 


1 d^S 
-f - 


2 9^2 


{V-Vf-pV. (A9) 


In order for the expansion to be about the maximum of 
the distribution, we must require the linear order term 
to vanish, i.e., we have P{V) = /3, meaning equivalence 
of the microcanonical and the canonical compactivity de¬ 
finitions. Identifying the inverse of —d^SjdV^ with the 
variance, our distribution takes the form 


P{V,P,N) = 


^S{V,N)-pV-{y-vY/2a^ 


(AIO) 


Evaluating this at two different compactivities and taking 
the logarithm of the ratio, we find (denoting the mean 
volumes by Vi and V 2 again) 


Q(V) = In 


P(V,P,,N) (E-Ei)2 , (E-E 2 )' 


P{V,P2,N) 2aj 

- PiVi + S{VuN) + P 2 V 2 - S{V2, N) 

z{P2,n) 


2 a| 


In 


Z{Pi,N) ’ 


(All) 


which is nothing but ( |A1[ ) with an explicit expression for 
ln((T2/(Ti). But we have derived this as an approximation 
to the distribution (181 from which we obtain Eq. (20) 


for Q. If we substitute P 2 for /3o in that equation, we see 
that the following (non-trivial) approximation holds in 
sufficiently large systems (close to the “thermodynamic 
limit”), as long as the distributions overlap significantly 


(which of course becomes less likely with increasing sys¬ 
tem size): 

(Y _ Y )2 (V - Vo)^ 

- (Pi - P2)V « ^ 


2cr^ 


2ai 


- piVi + S(Vi,N) + P 2 V 2 - S{V2, N) . (A12) 

Hence, the sum on the right-hand side that is quadratic 
in E is a good approximation to the sum on the left-hand 
side that is linear in E. As we have shown by this small 
calculation, the overlapping-histogram method will give, 
for such a system, the correct linear dependence on E, 
despite the fact that the central part of the distribution 
is well approximated by a Gaussian. In simulations, this 
behavior might be distinguished from Gaussian distribu¬ 
tions not having the statistical mechanical origin postu¬ 
lated by the Edwards theory through verification that 
the tails of the simulated distributions, i.e., their behav¬ 
ior for E values, where the quadratic approximation (A9) 
breaks down, are not Gaussian. This was done for our 
simulations via the chi-square test mentioned in Sec. |Vj 


Appendix B: Kernel density estimation 

A method to determine a continuous probability den¬ 
sity from a data sample is the kernel density estimation 
method (KDE) [53l[54j. If xi,X 2 , ■■■, Xn are sampled data, 
the kernel density estimation of the probability density 
P(x) at the point x is defined as 


1 




2=1 


X — Xj 


(Bl) 


where K(x) is the kernel which must be a non negative 
function that satisfies 


/ OO 

Ax K{x) = 1. 

-00 


(B2) 


and h IS a. smoothing parameter called the bandwidth. 
Possible kernels are for example the normal kernel 


K(x) = — exp(-a;2/2) , 


the Gauchy kernel 


K = 


I 


7 r(I -I- x"^) ’ 

the Epanechnikov kernel 

K(t) = ifa:e[-I;I] 

10 elsewhere 

or even the rectangular kernel: 

"l if a: G [-1/2; 1/2] 


K{t) = 


0 elsewhere 


(B3) 


(B4) 


(B5) 


(B6) 
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The latter is equivalent to a histogram with bin width h. 
While it can be shown that the Epanechnikov kernel is 
optimal in the sense that it minimizes the mean squared 
error between the estimated and the real probability dis¬ 
tributions, we used a normal kernel since it allows to 
make a good estimation of the optimal bandwidth h. In 
general, the optimal bandwidth can only be calculated 
if one knows the underlying probability density, but this 
density is unknown. In practice, therefore, Silverman’s 
rule of thumb is commonly used. Under the assumption 
that the underlying probability distribution is Gaussian 
and if a Gaussian kernel is used, the optimal bandwidth 
is 


L06.n-‘/» . (B7) 


where a is the standard derivation of the sample. It turns 
out that this bandwidth is also a reasonable choice in 
practical situations, if the underlying distribution is not 
Gaussian. We note that for small data sets the choice of 
the kernel may have a significant influence on the qual¬ 
ity of the fit. However, if the data sample becomes big 
enough, all kernels leads to almost the same results ex¬ 
cept for the far tail of the distribution, where no data 
points are available. In this region, the kernel itself spec¬ 
ifies the decay of the distribution. As it is shown exem- 
plarily in Figure the choice of the kernel is not crucial 
for our data samples. The results obtained with the op¬ 
timal Epanechnikov Kernel and the results achieved with 
the normal kernel are practically indistinguishable. The 
box kernel which is equivalent to a shifted histogram is a 
little bit more more irregular. We preferred the normal 
kernel, in order to avoid an ad hoc choice of the kernel 
bandwidth. 
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